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Abstract 

Motivated by recent progress on the interplay between graph the- 
ory, dynamics, and systems theory, we revisit the analysis of chemical 
reaction networks described by mass action kinetics. For reaction net- 
works possessing a thermodynamic equilibrium we derive a compact 
formulation exhibiting at the same time the structure of the complex 
graph and the stoichiometry of the network, and which admits a di- 
rect thermodynamical interpretation. This formulation allows us to 
easily characterize the set of equilibria and their stability properties. 
Furthermore, we develop a framework for interconnection of chemical 
reaction networks. Finally we discuss how the established framework 
leads to a new approach for model reduction. 

1 Introduction 

Large-scale chemical reaction networks arise abundantly in bio-engineering 
and systems biology. A very simple example involving only two chemical 
reactions in the glycolytic pathway is given by the following two coupled 
reactions 

Acetoacetyl ACP + NADPH + H + ^ D-3-Hydroxybutyryl ACP + NADP+ 
D-3-Hydroxybutyryl ACP ^ crotonyl ACP + H 2 0. (1) 
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The compounds Acetoacetyl ACP, NADPH, H+, D-3-Hydroxybutyryl ACP, 
NADP, crotonyl ACP and H2O involved in the reaction network are called 
the chemical species of the network. The most basic law prescribing the 
dynamics of the concentrations of the various species involved in a chemical 
reaction network is the law of mass action kinetics, leading to polynomial 
differential equations for the dynamics of each species. Large-scale networks 
thus lead to a high-dimensional set of coupled polynomial differential equa- 
tions, which are usually difficult to analyze. 

In order to gain insight into the dynamical properties of complex re- 
action networks it is important to identify their underlying mathematical 
structure, and to express their dynamics in the most compact way. In line 
with the recent surge of interest in network dynamics at least two aspects 
should be fundamental in such a mathematical formulation: (1) a graph 
representation, and (2) a specific form of the differential equations. 

The graph representation of chemical reaction networks is not immediate, 
since a chemical reaction (the obvious candidate for identification with the 
edges of a graph) generally involves more than two chemical species (the 
most basic candidate for identification with the vertices). We will follow 
an approach that has been initiated and developed in the work of Horn & 
Jackson [HI [13] and Feinberg, see e.g. [HUH], by associating the complexes 
of the chemical reaction network (i.e., the left- and right-hand sides of the 
reactions) with the vertices of the graplfl- The resulting directed graph, 
called the complex graph in this paper, is characterized by its incidence 
matrix. The expression of the complexes in the chemical species defines an 
extended stoichiometric matrix, called the complex stoichiometric matrix, 
which is immediately related to the standard stoichiometric matrix through 
the incidence matrix of the complex graph; see e.g. [21], [2]. 

In order to derive a specific form of differential equations we will start 
with the basic assumption that the chemical reaction rates are governed by 
mass action kinetics. As an initial step we then derive a compact form of 
the dynamics involving a non-symmetric weighted Laplacian matrix of the 
complex graph. The basic form of these equations can be already found 
in the innovative paper by Sontag [23] . The main part of the paper is 
however devoted to a subclass of mass action kinetics chemical reaction net- 
works, where we assume the existence of a thermodynamical equilibrium, 
or equivalently, where the detailed balance equations are assumed to admit 
a solution. We will call such chemical reaction networks balanced chemi- 



1 For alternative approaches based on species-reaction graphs and Petri-nets we refer 
to e.g. [61 0E]. 



cal reaction networks. Balanced chemical reaction networks are necessarily 
reversible but involve additional conditions on the forward and reverse re- 
action rate constants (usually referred to as the Wegscheider conditions; see 
|12| ) . For such balanced chemical reaction networks we will be able to derive 
a particularly elegant form of the dynamics, involving a symmetric weighted 
Laplacian matrix of the complex graph. Furthermore, it turns out that this 
form has a direct thermodynamical interpretation, and in fact can be re- 
garded as a graph-theoretic version of the formulation derived in the work 
of Katchalsky, Oster & Perelson [T9lfT8]. 

The obtained form of the equations of balanced chemical reaction net- 
works will be used to give, in a very simple and insightful way, a characteri- 
zation of the set of equilibria, and a proof of the asymptotic convergence to 
a unique thermodynamic equilibrium corresponding to the initial condition 
of the system. Similar results for a different class of mass action kinet- 
ics reaction networks, in particular weakly reversible networks with zero- 
deficiency, have been derived in the fundamental work of Horn |13) and 
Feinberg [10\ [TT] , which was an indispensable source of concepts and tools 
for the work reported in this paper. 

Subsequently we show how this form of the dynamics of chemical re- 
action networks can be extended to open reaction networks; i.e., networks 
involving an influx or efflux of some of the chemical species, called the bound- 
ary chemical species. Furthermore, we show how corresponding outputs can 
be defined as the chemical potentials of these boundary chemical species 
and how this leads to a physical theory of interconnection of open reaction 
networks, continuing upon the work of Oster &: Perelson |18| 120], 

In the final part of the paper we make some initial steps in showing 
how the derived form of balanced reaction networks can be utilized to de- 
rive in a systematic way reduced-order models. This reduction procedure is 
based on the Laplacian matrix of the complex graph, and draws inspiration 
from a similar technique in electrical circuit theory, sometimes referred to as 
Kron reduction. The model reduction procedure leads to an approximating 
reduced-order model which is again a mass action kinetics balanced chemical 
reaction network. 

Interestingly enough, the derived form of the equations of chemical reac- 
tion networks involving a weighted Laplacian matrix of the complex graph 
is reminiscent of the dynamics of the standard consensus algorithm of multi- 
agent systems, see e.g. [IT] . In fact, the form of the stability analysis in both 
cases is quite similar. Furthermore, the relation between the initial concen- 
tration of chemical species and the unique final equilibrium concentration is 
closely related to the recent concept of x-consensus in [5]. 



The paper is organized as follows. In Section [21 we summarize the 
mathematical structure of a network of chemical reactions described in e.g. 
[HI [131 EH EU EI] ; see also [2] for a clear account. In Section O we recall the 
law of mass action kinetics, and derive a general graph-theoretic formula- 
tion which is basically already contained in the innovative paper by Sontag 
[23]. Next, for the important subclass of chemical reaction networks admit- 
ting a thermodynamic equilibrium (balanced reaction networks) we derive a 
particularly simple formulation involving a symmetric weighted Laplacian 
matrix for the complex graph. Before providing the thermodynamical in- 
terpretation of this formulation, we discuss the graph-theoretic treatment 
of the complex graph, largely following the exposition in [21], and its im- 
plications for the established formulation of balanced reaction networks. In 
Section [4] we utilize the developed formulation to derive a simple character- 
ization of the set of equilibria, and we give a Lyapunov analysis to show the 
convergence to a unique equilibrium depending on the initial concentration 
vector. At the end of this section we comment on the relation with consen- 
sus algorithms. In Section [5] the framework is extended to reaction networks 
with inflows and outflows, and we discuss interconnection of such reaction 
networks through shared boundary chemical species. The initial steps in the 
application of the established formulation of balanced reaction networks to 
model reduction are discussed in Section [6l Finally, Section [7] contains the 
conclusions, and some topics of ongoing research. 

Notation: The space of n dimensional real vectors is denoted by M. n , and 
the space ofmxn real matrices by W nxn . The space of n dimensional real 
vectors consisting of all strictly positive entries is denoted by K™ and the 
space of n dimensional real vectors consisting of all nonnegative entries is 
denoted by R". The rank of a real matrix A is denoted by rank A. Given 
a\,. . . ,a n 6l, diag(ai, . . . , a n ) denotes the diagonal matrix with diagonal 
entries oi, . . . , a n ; this notation is extended to the block-diagonal case when 
a±, . . . ,a n are real square matrices. Furthermore, ker A and span A denote 
the kernel and span respectively of a real matrix A. If U denotes a linear 
subspace of M m , then U 1 - denotes its orthogonal subspace (with respect to 
the standard Euclidian inner product). l m denotes a vector of dimension m 
with all entries equal to 1. The time-derivative ^jjr(t) of a vector x depending 
on time t will be usually denoted by x. 

Define the mapping Ln : — > M. m , x i-4 Ln(x), as the mapping whose 
i-th component is given as (Ln(x)) i := ln(xj). Similarly, define the mapping 
Exp : lR m —7- R™, x i — y Exp(x), as the mapping whose i-th component 
is given as (Exp(x)) i := exp(xj). Also, define for any vectors x,z € M. m 



the vector x ■ z € M m as the element-wise product (x • z) i := XiZi, i = 

I, 2, ... ,m, and the vector | e M m as the element-wise quotient (|). := 

I I , i = 1, ■ ■ ■ ,m. Note that with these notations Exp(x + z) = Exp(x) • 
Exp(z) and Ln(x • z) = Ln(x) + Ln(z), Ln (|) = Ln(x) — Ln(z). 



2 Chemical reaction network structure 

In this section we will survey the basic topological structure of chemical re- 
action networks. First step is the stoichiometry expressing the conservation 
laws of chemical reactions. A next innovative step was taken in the work of 
Horn &: Jackson and Feinberg |14| [T3l [TU1 [TT] by defining the complexes of 
a reaction to be the vertices of a graph. We will summarize these achieve- 
ments in a slightly more abstract manner, also making use of the exposition 
given in [21]; see also [2] for a nice account. 



2.1 Stoichiometry 

Consider a chemical reaction network involving m chemical species (metabo- 
lites), among which r chemical reactions take place. The basic structure 
underlying the dynamics of the vector x 6 M.™ of concentrations Xi , i = 
1, • • • , m, of the chemical species is given by the balance laws 

x = Sv, (2) 

where S is an m x r matrix, called the stoichiometric matrix. The elements 
of the vector »£l r are commonly called the (reaction) fluxes. The stoichio- 
metric matrix S, which consists of (positive and negative) integer elements, 
captures the basic conservation laws of the reactions. For example, the 
stoichiometric matrix of the two coupled reactions involving the chemical 
species X\ , X2 , X3 given as 

X x + 2X 2 ^ X 3 ^ 2Xi + X 2 



is 



S 



-1 2 
-2 1 
1 -1 



while the first reversible reaction in (TjQ) 
Acetoacetyl ACP + NADPH + H+ ^ D-3-Hydroxybutyryl ACP+NADP" 



involving the species {Acetoacetyl ACP, NADPH, H + , D-3-Hydroxybutyryl 
ACP, NADP + , crotonyl ACP, H2O} has the stoichiometric matrix 



S=[-l 



-1 



-1110 



T 



In many cases of interest, especially in biochemical reaction networks, chem- 
ical reaction networks are intrinsically open, in the sense that there is a 
continuous exchange with the environment (in particular, flow of chemical 
species and connection to other reaction networks). This will be modeled by 
a vector Vb € M fe consisting of b boundary fluxes Vb, leading to an extended 
model 



x = Sv + S{,Vb 



(3) 



Here the matrix consists of mutually different columns whose elements 
are all except for one element equal to 1 or —1. The non-zero elements 
correspond to boundary fluxes for part of the chemical species, called the 
boundary chemical species, with +1 denoting an uptake flux and —la de- 
mand flux. 

Example 2.1. The open reaction network 

Ax + 2X 2 + 2X3 ^ 2X 4 + 2X 5 + 2X 6 ' 

- *3 > (4) 

A 6 ^ J 

has the extended stoichiometric matrix 



"-1 
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-2 








-2 


1 





2 








2 








_ 2 





-1 



[S Sb 



with boundary chemical species X3 and Xq . 

Although the stoichiometry does not fully determine the dynamics of 
the chemical reaction network (for this to be the case the vector of fluxes v 
needs to be expressed as a reaction rate v = v(x) ) it already contains useful 
information about the network dynamics, independent of the precise form 
of the reaction rate v(x). For example, if there exists an m-dimensional 
row-vector k such that 



kS = 



then the quantity kx is a conserved quantity for the dynamics x = Sv(x) for 
all possible reaction rates v = v(x). Indeed, ^kx = kSv(x) = 0. If k G M.™, 
then the quantity kx is commonly called a conserved moiety. Geometrically, 
for all possible fluxes the solutions of the differential equations x = Sv{x) 
starting from an initial state xo will always remain within the affine space 
xq + imS. In case of an open chemical reaction network j^kx = modifies 
into 

^-kx = kS b v b , 
at 

expressing that the time-evolution of the quantity kx only depends on the 
boundary fluxes v b . 

2.2 The complex graph 

The network structure of a chemical reaction network cannot be directly 
captured by a graph involving the chemical species (since generally there 
are more than two species involved in a reaction). Instead, we will follow an 
approach originating in the work of Horn & Jackson [141 US] an d Feinberg 
[101 Hlj - introducing the space of complexes. The set of complexes of a 
chemical reaction network is simply defined as the union of all the different 
left- and right-hand sides (substrates and products) of the reactions in the 
network. 

For example, the reaction network in (pQ) entails four complexes, namely 
the substrates Acetoacetyl ACP+ NADPH +H + and D-3-Hydroxybutyryl 
ACP+NADP + , and the products D-3-Hydroxybutyryl ACP and crotonyl 
ACP +H2O. In general, a product complex of one reaction may be the 
substrate complex of another. Thus, the complexes corresponding to the 
two reactions 

2Xi ^ Xy + 2X 2 ^±X 2 +X 3 

are 2Xi, Xi+A"2 and X2+X3. A complex may also be the product/substrate 
complex of more than one reaction. 

The expression of the complexes in terms of the chemical species is for- 
malized by an m x c matrix Z, whose p-th column captures the expression 
of the p-th complex in the m chemical species. For example, the column 
expressing the complex X\ + 2X2 + 2X3 in the chemical reactions Q is 



[1 2 2 0] 



Note that by definition all elements of the matrix Z are non-negative inte- 
gers. 

Since the complexes are left- and right-hand sides of the reactions they 
can be naturally associated with the vertices of a directed graph, with edges 
corresponding to the reactions. The complexes on the left-hand side of 
the reactions are called the substrate complexes and those on the right-hand 
side of the reactions are called the product complexes. Formally, the reaction 
a ^ 7r between the cr-th and the 7r-th complex defines a directed edge with 
tail vertex being the <r-th complex and head vertex being the 7r-th complex. 
The resulting graph is called the complex graph. 

Recall, see e.g. [3], that any directed graph is defined by its incidence 
matrix B. This is an c x r matrix, c being the number of vertices and r 
being the number of edges, with (p, j)-th element equal to —1 if vertex p is 
the tail vertex of edge j and 1 if vertex p is the head vertex of edge j, while 
otherwise. Obviously there is a close relation between the matrix Z and 
the stoichiometric matrix S, which is expressed as 



with B the incidence matrix of the complex graph. For this reason we will 
call Z the complex stoichiometric matrix. Hence the relation x = Sv between 
the fluxes v through the chemical reaction network and the time-derivative 
of vector of chemical species x, cf. ([2]), can be also written as 



with the vector Bv belonging to the space of complexes R c . 

3 The dynamics of mass action kinetics chemical 
reaction networks 

In this section we will derive a compact form for the dynamics of a chemical 
reaction network, whose reactions are described by mass action kinetics. 
After deriving a general form in Section 3.1 we will focus in the rest of the 
section on balanced chemical reaction networks. The resulting form of the 
dynamical equations for balanced reaction networks will be fundamental to 
the analysis of this dynamics in the subsequent sections. 

3.1 The general form of mass action kinetics 

The dynamics of the concentration vector x (or equivalently, in case of a 
fixed volume, the vector n of mole numbers) is given once the internal fluxes 



S = ZB, 



(5) 



x = ZBv 



(6) 



v are specified as a function v = v(x) of x, defining the reaction rates. The 
most basic model for specifying the reaction rates is mass action kinetics, 
defined as follows. Consider the single reaction 

X\ + X2 ^ X3 , 

involving the three chemical species X\, X2, X% with concentrations x\,X2,x$. 
In mass action kinetics the reaction is considered to be a combination of the 
forward reaction X\ + X2 — 1 X3 with forward rate equation v ozw (xi, X2) = 
k iorw xiX2 and the reverse reaction X\ + X2 t— X3, with rate equation 
v rev (x^) = k rcv X3. The constants /c forw , k rcv are called respectively the 
forward and reverse reaction constants. The net reaction rate given by mass 
action kinetics is thus 

V(X U X2,X 3 ) = V*°™(X 1 ,X2) ~ V" CV {X 3 ) = k l ° rW Xl X2 ~ ^3 

In general, the mass action reaction rate of the j-th reaction of a chemical 
reaction network, from a substrate complex Sj to a product complex Vj, is 
given as 

m m 

vj & = k?™ n x t ^ ~ ^ n x, ^ , (?) 

i=l i=l 

where Zj p is the (i,p)-th. element of the complex stoichiometric matrix Z, 
and kj W ,kj SV > are the forward/reverse reaction constants of the j'-th 
reaction, respectively. Without loss of generality we will throughout assume 
that for every j the constants £^ orw , /c^ cv are not both equal to zero (since in 
this case the j-th reaction is not active). 

This can be rewritten in the following way. Let Zg i and Z-p. denote 
the columns of the complex stoichiometry matrix Z corresponding to the 
substrate and the product complexes of the j-th reaction. Using the mapping 
Ln : IRj_ — > M c as defined at the end of the Introduction, the mass action 
reaction equation (J7|) for the j-th reaction from substrate complex Sj to 
product complex Vj can be rewritten as 

Vj (x) = kf TW exp (Zj. Ln(x)) - kf v exp (z£ Ln(x)) . (8) 



Based on the formulation in (|8|), we can describe the complete reaction 
network dynamics as follows. Let the mass action rate for the complete set 

rp 

of reactions be given by the vector v(x) = [v±(x) ■ ■ ■ v r (x)~\ . For every 
a, 7T G {!,■■■ ,c}, denote by a a7r = kf\ = kf™ if (a,n) = (Sj,Vj), 



j € {1, . . . , r} and a a7T = elsewhere. Define the weighted adjacency matrix 
A of the complex graprj^l as the matrix with (a, 7r)-th element a„, where 
a, ir G {1, • • • , c}. Furthermore, define the weighted Laplacian matrix L as 
the c x c matrix 



L := A - A (9) 

where A is the diagonal matrix whose (p, p)-th element is equal to the sum 
of the elements of the p-th column of A. Then it can be verified that the 
vector Bv{x) for the mass action reaction rate vector v(x) is equal to 

Bv{x) = -LExp (Z T Ln(x)) , 

where the mapping Exp : R c — > R+ has been defined at the end of the 
Introduction. Hence the dynamics can be compactly written as 

x = -ZLExp (Z T Ln(x)) (10) 

A similar expression of the dynamics corresponding to mass action kinetics, 
in less explicit form, was already obtained in [23]. 



3.2 Balanced mass action kinetics 

In the rest of this paper we will focus on an important subclass of mass 
action chemical reaction networks. We need the following basic notions. 

Definition 3.1. A vector of concentrations x* E R™ is called an equilib- 
rium^ for the dynamics x = Sv(x) if 

Sv(x*) = 

Furthermore, x* E is called a thermodynamic equilibrium if 
v{x*) = 

Clearly, any thermodynamic equilibrium is an equilibrium, but not nec- 
essarily the other way around (since in general S = ZB is not injective). 

2 Strictly speaking, A is the adjacency matrix not of the complex graph, but of a directed 
augmented graph. Each j-th edge of the complex graph is replaced by two directed edges: 
one corresponding to the forward reaction with weight kj° rw and one corresponding to the 
reverse reaction with weight . 

' Sometimes called a kinetic equilibrium, in order to distinguish it from a thermody- 
namical equilibrium. 



3.2.1 The existence of thermodynamic equilibria 

Necessary and sufficient conditions for the existence of a thermodynamic 
equilibrium can be derived in the following linear-algebraic way following 
|12j . These conditions are usually referred to as the Wegscheider conditions, 
generalizing the classical results of |24j . 

Consider the j-th reaction from substrate Sj to product Vj, described 
by the mass action rate equation (see ©) 



Vj (x) = kf™exp (Zj.Ln(x)) - fcf v exp (2$.Ln(x)) 

Then x* € W£ is a thermodynamic equilibrium, i.e., v(x*) = 0, if and only 
if 

fcf™exp(zJ.Ln(z*)) =^ ev exp(Z^Ln(x*)), j = 1, • • • ,r (11) 
The equations (jlip . sometimes referred to as the detailed balance equations, 



7-eq 

reaction as (assuming fc^ ev ^ 0) 



can be rewritten as follows. Define the equilibrium constant K- of the j-th 



3 

r,forw 

K eq := J— (12) 

K j 



Then the detailed balance equations (jlip are seen to be equivalent to 

Kf = exp (z£.Ln(x*) - Zj.Ln(x*)) , j = 1, • • • , r (13) 

Collecting all reactions, and making use of the incidence matrix B of the 
complex graph, this amounts to the vector equation 

K eq = Exp (B T Z T Ln(x*)) = Exp (,S T Ln(x*)) , (14) 

where K eq is the r-dimensional vector with j-th element K eq ,j = 1, • • • , r. 

From here it is easy to characterize the existence of a thermodynamic 
equilibrium. 

Proposition 3.2. There exists a thermodynamic equilibrium x* S WV: if 
and only if h!° TW > 0, kj ev > 0, for all j = 1, • • • , r, and furthermore 

LnCFT 9 ) €imS T (15) 

Proof. Clearly k { ° TW > 0, kj cv > 0, j = 1, ■ ■ ■ , r, is a necessary condition for 
the existence of a thermodynamic equilibrium. The existence of a vector 
Ln(x*),x* £ R+, satisfying (fill is obviously equivalent to (fT5l) . ■ 



It also follows that once a thermodynamic equilibrium x* is given, the 
set of all thermodynamic equilibria is described as follows. 



Proposition 3.3. Let x* € be a thermodynamic equilibrium, then the 
set of all thermodynamic equilibria is given by 

£ := {x** 6 [ S T Ln (x**) = S T Ln (x*)} (16) 

For ease of exposition we will henceforth refer to mass action chemical 
networks possessing a thermodynamic equilibrium x* € WV as balanced mass 
action chemical networks: 

Definition 3.4. A chemical reaction network x = Sv(x) governed by mass 
action kinetics is called balanced if there exists a thermodynamic equilibrium 
x* e R+. 

As stated in Proposition 13. 2[ a necessary condition for the existence of 
a thermodynamic equilibrium is the fact that for all reactions the forward 
and reverse reaction constants are both strictly positive. Thus all reactions 
of a balanced reaction network need to be at least reversible. 

Remark 3.5. Usually the conditions \15\) are rewritten in the following 
more constructive form. By basic linear algebra A15\) is satisfied if and only 
if for all row-vectors a satisfying o~S T = we have 

J>,lni<7 = (17) 
i=i 

Putting back in the definition of the equilibrium constants U2\) this is seen 
to be equivalent to the usual form of the Wegscheider condition^ 

n(fcf™) <Tj =f[(k™p. (is) 

i=i j=i 

4 Note that (|18|) only needs to be checked for a maximally independent set of row- 
vectors a satisfying aS T = 0. Furthermore, by writing this as aS T = aB T Z T = this 
can be related to the topological structure of the complex graph and to the deficiency of 
the network; cf. [12]. Furthermore, the Wegscheider conditions admit a thermodynamical 
interpretation, since they are intimately related to the passivity properties of the network 
and the independence of the increase of Gibbs' energy of the path in the reaction network; 
see also [181119] . 



3.2.2 The standard form of balanced mass action reaction net- 
works 



For balanced mass action chemical reaction networks we can further rewrite 
the dynamics (|1U|) in the following useful way. Let x* G be a thermody- 
namic equilibrium, i.e., v(x*) = 0. Consider the rewritten form (|13p of the 
'detailed balance' equations. These equations allow us to define the balanced 
reaction constant of the j-th reaction as 

Kj(x*) : = £;j orw exp (zJXn(x*)) = kf v exp (z£Xn(x*)) , j = 1, ••• ,r 

(19) 

Then the mass action reaction rate of the j-th. reaction can be written as 

r T - ( x 



Vj (x) = Kj (x*) [exp \Z Sj Ln j j - exp (Zp.Ln 

where for any vectors x,z € R m the quotient vector | € R m is defined 
elementwise (see the end of the Introduction). Defining the r x r diagonal 
matrix of balanced reaction constants as 

tC(x*) := diag(/ci(ar*), • • • , K r (x*)) (20) 

it follows that the mass action reaction rate vector of a balanced reaction 
network can be written as 

v ( x ) = -/C(x*)S T Exp (z T Ln ( — 

\ \x* 

and thus the dynamics of a balanced reaction network takes the form 

x = -ZBJC(x*)B T Exp(z T Ln(^J^ , K(x*) > (21) 

This form will be the starting point for the analysis of balanced chemical 
reaction networks in the rest of this paper. 

The matrix BK.(x*)B T in $FQ) again defines a weighted Laplacian matrix 
for the complex graph, with weights given by the balanced reaction constants 
K%(x*), ■ ■ ■ ,K r (x*). Note that this is in general a different weighted Laplacian 
matrix than the one obtained before, cf. ©. In particular, a main difference 
is that the weighted Laplacian BK{x*)B T is necessarily symmetric. Among 



5 Although these constants appear in the literature, see e.g. [18(119). we could not find 
a standard name for them. Note that Kj(x*) > is defined as the common value of the 
forward and reverse reaction rate of the j-th reaction at thermodynamic equilibrium. 



others, cf. [3], this implies that the Laplacian BK{x*)B T is in fact inde- 
pendent of the orientation of the graph. Thus we may replace any reaction 
S f± V by V f± S without altering the Laplacian BK.(x*)B T , in agreement 
with the usual understanding of a reversible reaction network. The sym- 
metrization of the Laplacian has been accomplished by the modification of 
Ln(x) into Ln (:p-), and using the assumption that x* is a thermodynamic 
equilibrium. 

Note that JC(x*), and therefore the Laplacian matrix BK(x*)B T , is in 
principle dependent on the choice of the thermodynamic equilibrium x*. In 
the next Section \3. 31 we will see that actually this dependence is minor: for a 
connected complex graph the matrix IC(x*) is, up to a positive multiplicative 
factor, independent of the choice of the thermodynamic equilibrium x*. 



3.3 The linkage classes of the complex graph 

The complex graph provides a number of tools for the analysis of reaction 
networks. Recall that for any directed graph [3] 

rank B = rank L = c - £, (22) 

where c is the number of vertices of the graph, and I is equal to the number 
of component^ of the complex graph , the linkage classes in the terminology 
of |14 |. 110 ^ fTT] . Furthermore, if there is one linkage class in the network (i.e., 
the graph is connected and ranki? = rankL = c — 1), then 

ker L = ker B T = span l c , 

where as before l c is the c-dimensional vector with elements all equal to 1. 

In general, if the reaction network has i linkage classes then the network 
can be decomposed as follows. Let the p-th linkage class have r p reactions 
between c p complexes. Partition Z, B and fC matrices according to the 



6 A directed graph is connected if there is a path (a number of un-oriented edges) 
between every two distinct vertices of the graph. The components of a directed graph are 
the maximal connected subgraphs. 



various linkage classes as 



Z = 



[Zi z 2 



Zt] 



B x 
B 2 













B = 



(23) 












Bi-\ 






Be 



K{x*) = diag(/Ci(x*),/C 2 (x*),...,/Q(x*)) 
where for p = 1, ...,£, Z p G M™ XCp , 5 p G M c P xr p and /C p G K+ pXrp denote, 



respectively, the complex stoichiometric matrix, incidence matrix and the 
diagonal matrix of balanced reaction constants for the p-th linkage class. 
Furthermore, S p = Z p B p is the stoichiometric matrix of the p-th. linkage 
class. It follows that equation (pZTj) can be expanded as 



expressing the contributions of each linkage class to the chemical reactionaj. 

The above expressions also yield the following alternative characteri- 
zation of the set of thermodynamic equilibria £ given in (|16p . Since for a 
connected graph ker B T = span l c , the set £ can be equivalently represented 
as, writing S T = B T Z T , 



For non-connected graphs the same formula holds for every component. 
This observation has the following useful consequence. 

Proposition 3.6. Consider the definition of the diagonal matrix )C(x*) of 
balanced reaction constants in 120\) . and the formulation of the dynamics 
given in h21\) for a particular thermodynamic equilibrium x* . Assume that 

7 Note that this yields a conserved moiety for the reaction network. In fact, the law 
of conservation of mass states that for any complex stoichiometric matrix Z there exists 
u € R™, such that Zju G spanl Cp for p = 1, . . . ,£. Since Bjl Cp = 0, this implies that 
u T x is a conserved moiety for x = ZBv{x), for all forms of the reaction rate v(x). 




(24) 



P =i 



£ = {x** G M+ | Z T Ln(x**) - Z T Ln(x*) G spanlj 



(25) 



the complex graph is connected. Then for any other thermodynamic equilib- 
rium x** there exists a positive constant d** such that 

JC(s~) = d**K(x*), Exp(Z r Ln (JL) ) = ^Exp(Z r Ln (JL) ) 

Furthermore, for a non-connected complex graph there exist positive con- 
stants d**,p = 1, • • • ,£, such that for each p-th linkage class the diagonal 
matrix K, p (x*) given in $23\) satisfies 

K p {x**) = d* p *K, p {x*\ P = i,---,£ 

Proof. Assume that the complex graph is connected. Then by (|25p there 
exists a constant c** such that Z T Ln(x**) = Z T hn(x*) + c**l. Hence by 
the definition of the balanced reaction constants ka in ([190 we have 

Kj(x**) = kf TW exp (z^.Ln(x**)\ = kf IW exp (z%.Ln(x*) + c**) 
= d**kf rw exp (z%.Ln(x*)\ =d**Kj(x*), 

with d** := exp(c**) > 0, for each j = 1, • • • ,r. The rest follows easily. ■ 

Hence the weighted Laplacian matrix BK,(x*)B T is, up to a multiplica- 
tive factor, independent of the choice of the thermodynamic equilibrium x*. 
The properties of this matrix (e.g., its eigenvalue distributiord) therefore 
may serve as an important indicator of the network dynamics. In Section [6] 
we will exploit the Laplacian matrix for model reduction purposes. 

3.4 Deficiency 

An important notion to relate the structure of the complex graph to the 
stoichiometry, as introduced in the work of Feinberg [10] , is the notion of 
deficiency. 

Definition 3.7. The deficiency of a chemical reaction network with complex 
stoichiometric matrix Z and incidence matrix B is defined as 

5 := rank B — rank ZB = rank B — rank S > (26) 

A reaction network is said to have zero- deficiency if 5 = 0. 



8 Note that the set of eigenvalues of B)C(x*)B T consists of (with multiplicity equal 
to the number of components) and strictly positive real numbers, cf. [3]. 



Note that zero-deficiency is equivalent to 



kerZ n imB = 0, 

or to the mapping Z : imB C M c — > R m being injective. Hence in the 
zero-deficiency case there is a one-to-one correspondence between the rate 
vector i £ im5 C M m of chemical species x E WT: and the rate vector 
G imB C l c of complexes y € R c . Many chemical reaction networks are 
zero-deficient, although with growing complexity (especially in biochemi- 
cal networks), deficiency greater than zero is likely to occur. This is also 
expressed in the following proposition showing that for a reaction network 
with zero-deficiency all its linkage classes also have zero-deficiency, but not 
necessarily the other way around. 

Proposition 3.8. 121)1 Consider a chemical reaction network with £ linkage 
classes as in \2J$ . If the network has zero- deficiency then necessarily the 
linkage classes also have zero- deficiency. However, zero- deficiency of the 
linkage classes does not generally imply zero- deficiency of the total network. 
This does hold if and only if additionally 

i 

C\ im Z p B p = 
P =i 

where Z p and B p denote the columns of Z, respectively of B, in \2J$ corre- 
sponding to Z p , respectively B p . 



3.5 Thermodynamical viewpoint 

In this section we will indicate the thermodynamical interpretation of the 
quantities introduced before, and show how this suggests a Lyapunov func- 
tion which will be used in the next section. For more details regarding the 
thermodynamical approach to chemical reaction kinetics we refer to [19} [T8] . 

Recall that for an ideal dilute solution the chemical potential [Xi of chem- 
ical species i with mole number m is given by 

qn ■ 

Hi(xi) = tf + RTln(^) = fi° + RTln(xi) (27) 

with n° a reference potential, R a constant, T the temperature, V the vol- 
ume, and Xi = y the concentration. Equivalently, we have the inverse 
relation 

Xi = exp 



RT 



The m-dimensional vector \i with components fii is called the chemical po- 
tential vector, while the vector [x° with components ji° is called the reference 
chemical potential vector. 

Starting instead from the formulation of a balanced chemical reaction 
network in (|2ip . corresponding to the thermodynamic equilibrium x*, we 
may define the chemical potential vector \x and the reference chemical po- 
tential vector as 

H(x) := RTLn (Jj) , fj,° := -RTLn{x*) (28) 

We conclude that Ln(pr) = Ln(x) — Ln(x*) is, up to the constant RT, equal 
to the chemical potential vector, while — Ln(x*) is, up to this same constant, 
equal to the reference chemical potential vector. Recalling from (|14p the 
expression of the vector of equilibrium constants as K eq = Exp (S' T Ln(x*)), 
it follows that 

RTLn (K eq ) = (RT)S T Ln(x*) = -S T n° (29) 

An important role in equation (I21D is formed by the quantity j(x) := 
Z T /j,(x) = Z T Ln (-tt) , which we will call the complex thermodynamical 
affinity. Correspondingly, we refer to 7 := Z T fi° as the reference complex 
thermodynamical affinitvQ. The dynamics of a balanced reaction network 
is determined by this complex thermodynamical affinity 7, which acts as a 
'driving force' for the reactions. 

Remark 3.9. Consider on the other hand the vector a := S T n, known as 
( minv^\) the vector of thermodynamical affinities. Then for mass action 
kinetics It is not possible, see e.g. /TP), [iffi . to express the vector of reac- 
tion rates v(x) as a function of a (contrary to j). This is illustrated in 
JW$ by considering the reaction X\ ^ X2 with reaction rate v(xi,X2) = 
kf° rw xi — k rev X2- When the concentrations x\ and X2 are doubled, then so 
is the reaction rate v{x\,X2). However, the thermodynamical affinity —a 
remains the same. 

In thermodynamics the vector of chemical potentials is derived as the 
vector of partial derivatives of the Gibbs' free energy. This suggests to define 
the Gibbs' free energy as 

G{x) = RTx T hn (^j + RT (x* - x) T t m , 

9 Note that (RT)Ln(K eq ) = (RT)S T Ln(x*) = -B T -y°, and thus the equilibrium con- 
stants correspond to differences of the reference complex thermodynamical affinities of the 
substrate and product complexes. 

10 To remain in line with the usual sign convention in thermodynamics. 



where l m denotes a vector of dimension m with all ones. Indeed, as is 
immediately checked 



§§(«)- 

In the rest of this paper we will assume that RT = 1, or equivalently, we 
will (re-)define 

:= Ln (JL) , 7 ( x ) := z T ^x) = Z T Ln (£) 

(30) 

G(x) := x T Ln + (x* - xf l m 



It follows that the equations of a balanced chemical reaction network ([2T 
can be equivalently written as 



7T dG , A , x dG 



x = —ZRK(x*)B Exp Z 1 —(x) , a(x) = —(x) (31) 

\ ox J ox 

In the next section, we will employ G(x) as a Lyapunov functioi^l for the 
chemical reaction network. In particular we will prove that ^j? < 0; showing 
that G is non-increasing along solution trajectories. 

A geometrical interpretation of Equations (12 ip and (|3ip is as follows. 
Denote the dual space of the space of concentrations of chemical species 
M := M. m by M*. Similarly, denote the dual space of C := R c by C*, and the 
dual of the space of reaction rates 1Z = W by 1Z* . Define v* := B T 'Exp( , y) 
and y := Bv(x). All ingredients of the equation (|3ip are then summarized 
in the following diagram 

v G TZ — > yeC x£M 



K{x*) | | G{x) 

(32) 

v* € TZ* < — 7 G C* < — u£M* 
B T Z T 

o 

Exp 

which expresses the duality relations between all the variables involved. The 
concentration vector x and its time-derivative x are elements of the linear 

11 Note that the definition of G depends on the chosen thermodynamical equilibrium. 
Denoting the functions for different thermodynamic equilibria x* and x** by G* , respec- 
tively G", it is seen that G**(x) — G*(x) + x T [Ln(a;*) — Ln(x**)] + (x** — x*) T \ 7n . 



space M with as conjugate vector the chemical potential vector \i G A4* . 
They are related by the Gibbs' function G(x) as fi = QQ{x). Furthermore, 
the vector y is in the linear space C, with conjugate vector the complex 
affinity 7. The relations between y and x and dually between \x and 7 are 
given by x = Zy, respectively 7 = Z T 11. Also note that x = Zy = ZBv = 
Sv, where the vector of fluxes v is in the linear space 1Z, with conjugate 
vector v* := —UC(x*)) v € 1Z*. The added complication in the diagram 
is the map Exp : C* — > C* , which introduces a discrepancy between v* and 
a ■= -B T j = -S T fi. 

4 Equilibria and stability analysis of balanced re- 
action networks 

In this section we show that the set of equilibria of a balanced chemical 
reaction network are actually equal to the set of thermodynamic equilibria 
£, and we study their stability properties. 

4.1 Equilibria of balanced reaction networks 

Making use of the formulation of the dynamics of balanced reaction net- 
works in (|2ip . we give a simple proof of the statement that all equilibria 
of a balanced reaction network are actually thermodynamic equilibria, and 
thus given by (|16p . A similar result was obtained in the classical papers 
|13| , I14| , fTT] for a different class of chemical reaction networks (roughly speak- 
ing, weakly reversible networks of deficiency zero or deficiency one under 
additional conditions). 

Theorem 4.1. Consider a balanced chemical reaction network x = Sv = 
ZBv with m species and r reactions governed by mass action kinetics, with 
thermodynamic equilibrium x* , i.e., v(x*) = 0, described as in \21\) . Then 
the set of all equilibria is equal to the set £ := {x** £ W£ \ S T Ln(x**) = 
S T Ln(x*) = 0} of thermodynamic equilibria given in l\16)) . 

Proof. Denote for the j-th reaction the substrate complex by Sj and the 
product complex by Vj. Let as before Z$ 3 and Z-p j denote the columns of 
the complex stoichiometric matrix Z of the reaction network, corresponding 
respectively to substrate complex Sj and product complex Vj. Define as 
before 



H(x) = Ln , 7 (x) := Z T [i(x), 7^(2;) = Z^fi(x), 7^0*0 = Z^.fj,(x) 



Suppose x** is an equilibrium, i.e., ZBIC(x*)B T Exp (Z T Ln (%r)) = 0, or 
equivalently ZBJC(x*)B T Ex.p (Z T fi(x**)) = 0. Then also 

fi T (x**)ZBIC(x*)B T Exp (Z T fi(x**)) = 

Denoting the columns of B by b\ , • • • ,b r , we have 
r 

BK{x*)B T = ^K j {x*)b j b] 

3=1 

while 

= ^(x**) - Z S .) = 7 £.(***) - !%{x*% 

bjExp (Z T fi{x**)) = exp (t^.(x**)) - exp (75/2;**)) 

It follows that 

= n T (x**)ZB1C(x*)B T Exp (Z T n(x**)) 

= E -=i [ivM*) ~ 75,00] [exp (TPj (»**)) " exp ( 7 S, (***))] ^(x*) 

Since the exponential function is a strictly increasing function and nj(x*) > 
0, j = 1, • • • , r, this implies that all terms in the summation are zero and 
thus 

IVj ( x **) = 1S 3 (x**), j = 1, • • • , r 

exp (jv 3 (x**)) = exp (7^ (x**)) , j = 1, • • • , r 

which is equivalent to 

£ T 7 (x**) = B T Z T n{x**) = 0, B r Exp ( 7 (x**)) = 5 T Exp (Z T ^(x**)) = 

The first equality tells us that x** € £, and thus a thermodynamical equi- 
librium (as also follows from the second equality). ■ 

4.2 Asymptotic stability 

In the next theorem we show that G serves as a Lyapunov function for (|2ip , 
implying global asymptotic stability with respect to the set of equilibria £ 
given by (fl6l) . i.e., £ is globally attractive. 

Theorem 4.2. Consider a balanced mass action reaction network given by 
\21\) or, equivalently, by \31\) . Then for every initial condition x(0) E W?, 
the species concentration x converges for t —¥ 00 to £ . 



Proof. First we prove that the function G in ()30p satisfies 

G{x*) = 0, G(x)>0, Vx^x*, (33) 

and is proper, i.e., for every real c > the set {x 6 | G(x) < c} is 
compact. With regard to (|33|) we first note the following. Let Xj and x* 
denote the i-th elements of x and x* respectively. From the strict concavity 
of the logarithmic function 

z-l>\n(z), Vz £ i + , (34) 

x * 

with equality if and only if z = 1. Putting z = — in equation (1341) . we get 

x* - Xi+Xi In (j^) - °i 
with equality if and only if x, = x*. This implies that 

> 0, 



G(x) = x* - Xj + Xj In I — ^ j 



with equality if and only if = 1, • • • , to. Thus G has a strict 

minimum at x = x* and satisfies (|33p . Properness of G is readily checked. 
Next we will show that G(x) := ^g^-(x)Sv(x) = ^jf (x) satisfies 

G(x) < Vx G K™, (35) 

and 

G(x) = if and only if x e £ . (36) 

As in the proof of Theorem 14.11 let the j-th reaction be between substrate 
complex Sj and product complex Vj , and let Zs j and Zj>. denote the columns 
of the complex stoichiometric matrix Z corresponding to complexes Sj and 
Vj. Define as before 

H(x) = Ln (-j^J , 7 (x) := Z T fi(x), 7 5j (x) = Z^fi(x), yp s (x) = Z^.fj,(x) 

Then compute^! 

G{x) = ^?(x)x = -/x T (x)ZS/C(x*) J B T Exp(Z T /i(x)) 
= -7 T (x) J B/C(x*) J B T Exp( 7 (x)) 

= Ej=i bv 3 {x) - 75,0*0] [exp (TPj(aj)) - exp {is^x))] Kj{x*) 
< 0, 



See the proof of Theorem 14.11 for an analogous reasoning. 



(37) 



since Kj(x*) > for j = 1, ...,r, and the exponential function is strictly 
increasing. The "if" part in (|36p is trivial. For the "only if" part, the 
summand in the third line of (|37p is zero only if 75^ (x) — r ^-p i (x) = for 
every j. This is equivalent to having B T j(x) = 0. Thus, G(x) = only if 
B T j = B T Z T Ln (^r) = 0. It follows that G(x) = if and only if x € S. 

Suppose that the set MIp is not forward invariant with respect to (|3ip . 
and that Xi{t) = for some t and i £ {1, • • • ,m}. Without loss of gen- 
erality assume that the species with concentration Xi is taking part in at 
least one complex which is involved in a reaction. Let Ci denote the sub- 
set of complexes which contain Xi, i.e., Zid 7^ 0. Then it follows that 
Exp(Z^fi(x(t))) = and thus 

Xi {t) = -Z l BlC{x*)B T Exp(Z T 'fi(x)) > 0, 

where Z % is the i-th row vector of Z . This inequality is due to the fact 
that the terms corresponding to the positive i-th diagonal element of the 
weighted Laplacian matrix BK,{x*)B T are all zero, while there is at least 
one term corresponding to a non-zero, and therefore strictly negative, off- 
diagonal element of BK,{x*)B T . This implied! that must be forward 
invariant. 

Since G is proper (in M™) and the state trajectory x(-) remains in M™, 
([35]) implies that x(-) is bounded in Therefore, boundedness of x(-), 
together with equations ([35]) and (i36j) . imply that the species concentration 
x converges to £ by an application of LaSalle's invariance principle. ■ 

Remark 4.3. A similar reasoning for showing that G < was used in f^j, 
see also fT^[ [lu\ [77]/ : however making use of the convexity of the exponen- 
tial function instead of the weaker property that the exponential function is 
increasing as in our proof. In this respect, it can be noted that Theorem \4-%\ 
remains unaltered if we replace the dynamics 121\) by any equations of the 
form (not corresponding anymore to mass action kinetics) 

x = -ZBlC(x*)B T F(z T Ln(^)), 

with F : R c -> R c a mapping F(y 1 , ■■■ ,y c ) = diag(/i(yi), • • • , f c {y c )), where 
the functions /j, i = 1, • • • , c, are all monotonously increasing. 



A similar argument can be used directly for the mass action kinetics rate equations. 



4.3 Equilibrium concentration corresponding to an initial 
concentration 



In this section, we show that for every initial concentration vector xq G 
R!l, the solution trajectory of f)21 1) converges to a unique thermodynamic 
equilibrium in 6. Our reasoning is very similar to the proof of zero-deficiency 
theorem provided in [11] and is based on the following proposition in there. 
Recall from the Introduction that the product x • z G M m is defined as the 
element-wise product (x ■ z)i = X{Zi, i = 1, • • ■ , m. 

Proposition 4.4. Let U be a linear subspace o/R m , and let x*,xq G W?. 
Then there is a unique element \i G U ± , such that (x* ■ Exp(/x) — xq\ G U . 

Proof. See proof of [TTJ Proposition B.l, pp. 361-363]. ■ 

Theorem 4.5. Consider the balanced chemical reaction network \2 Then 
for every xq G there exists a unique x\ G £ such that the solution 

trajectory of \21\) starting from xq converges for t — > oo to x\. Hence there 
exists a surjective map \ '■ R+ — ^ £ that assigns to every initial state its 



asymptotic thermodynamic equilibrium \. 

Proof. With reference to Proposition 14.41 define U = span S, and observe 
that U 1 - = kerS T . Let x*,xq G M^. Then by Proposition 14.41 there exists 
a unique [i G ker S T such that x* ■ Exp(/i) — x$ G span 5. Define x\ := 
x* • Exp(//) G R!p. It follows that S T fi = S T Ln (ft) = 0, i.e., x x G £. Also, 
since x\ — xq G spanS", xi G {£ G M m | S,~ x o € span 5} which is an invariant 
set of the dynamics. Together with Theorem 14.2^ the state trajectory with 
initial condition xo converges to the unique equilibrium x\ G £. ■ 

4.4 Relation with consensus dynamics 

The form of the equations (|10p is reminiscent of consensus dynamics for 
multi-agent systems with non-symmetric communication. In fact, a stan- 
dard consensus algorithm in the vector y G M c can be written as [17] 

V — L cons y, 

where L C ons is the weighted Laplacian matrix of the communication graph. 
This weighted Laplacian is constructed similarly as in ([9]), with the differ- 
ence that A is now defined as the diagonal matrix whose (p, p)-th element 

14 The map \ i s an analytic function, following a similar argument as used in [231 The- 
orem 6] which deals with weakly-reversible zero-deficiency chemical networks. 




is equal to the sum of the elements of the p-th row of A, implying that the 
sum of the elements of every row of -L CO ns is zero. Note that in contrast for 
the weighted Laplacian L in ([9]) the sum of the elements of every column is 
zercf^l. Seen from this perspective, the equations (|2ip for balanced chemi- 
cal reaction networks, exhibiting the symmetric weighted Laplacian matrix 
BK,(x*)B T , correspond to the consensus algorithm for multi-agent systems 
with symmetric communication structure |17j . 

Finally, Theorem 14.51 is closely related to the property that the asymp- 
totic consensus value in consensus dynamics is a function of the initial state. 
In fact, the surjective map x '■ ^+ £ which assigns to every initial state 
its asymptotic equilibrium is similar to the x-function in the x- consensus 
algorithm of [5]. 

5 Chemical reaction networks with boundary fluxes 
and interconnection of reaction networks 

As discussed in Section 2.1, the interaction of a chemical reaction network 
with the environment or another chemical reaction network can be modeled 
by identifying a vector of boundary chemical species xj, € M 6 , which is a sub- 
vector of the vector x of concentrations of all the chemical species in the 
network. These boundary chemical species are the species that are subject 
to inflow or outflow boundary fluxes v b - The natural conjugate vector is (up 
to the constant RT) the vector of chemical potentials [i b = STfJ, E M. b of 
boundary potentials. Indeed, up to the factor RT, the product ^v b is equal 
to the power entering or leaving the chemical reaction network due to the 
influx or efflux of boundary chemical species. 

By combining (|3|) and (}2T|) this leads to the following equations 

x = -ZB1C(x*)B t Etcp (Z T Ln (4)) + S b v b , 

(38) 

which define an input-state-output system with inputfl v b and outputs fib- 
Let G be defined by equation ([30]) . Define as before 7(2;) := Z T Ln (^r). As 

15 Of course, another main difference between the two dynamics is the appearance of 
the complex stoichiometric matrix Z on the left of the right-hand side of the differential 
equation (| 10 j) . and the much more involved term Exp (Z T Ln(a;)) instead of simply y (see 
however [5]). 

16 Note however that the elements of Vb can be either inflows or outflows for the network. 



an immediate consequence we obtain the following energy balance 



d_ 

It 



G{x) = -^ T (x)B)C(x*)B T Exp( 7 (x)) + fj$v b 



(39) 



Furthermore, since the exponential function is increasing it is easily verified 
along the lines of the proof of Theorem 14.21 (see in particular (|37p ) . that 
7 T /CB T Exp(7) > 0, thus showing the passivity property 

jG < Ufa (40) 

5.1 Interconnection of chemical reaction networks 

The most basic way of interconnecting chemical reaction networks is through 
shared boundary chemical species. Consider for simplicity of notation the in- 
terconnection of two chemical reaction networks which have all their bound- 
ary chemical species in common; the generalization to the general case is 
straightforward. Let denote the incidence matrix of the complex graph, 
and Zi the complex stoichiometry matrix of network i = 1,2. The complex 
graph of the interconnected chemical reaction network is the direct union of 
the complex graphs of networks 1 and 2, with incidence matrix B given as 
the direct product 

\B X " 
B 2 



B 



(41) 



Note that the complex graph of such an interconnected network is not con- 
nected. In fact, if all the constituent networks are connected then their 
complex graphs are linkage classes of the interconnected network, cf. Sec- 
tion GTS 

The complex stoichiometry matrix of the interconnected network is more 
involved. Permute the chemical species x±,x 2 such that 

Zb2 



xi = (xi,x b i),x 2 = (x b2 ,x 2 ), Zi 



Zi 
Zti 



Zo 



z 2 



where Z b i,Z b2 are matrices with the same number of rows, equal to the 
number of shared boundary species 

Xb ■= x b i = x h2 € R b 

Then the stoichiometry matrix Z of the interconnected network is given as 

' ' Zi 

(42) 



Zbi 




Zb2 
Z2. 



Note that in general the property of zero- deficiency is not maintained under 
interconnection. In fact, this was already discussed in Proposition 13.81 in 
Section E3] following [21] . 

Remark 5.1. It may occur that one or more columns of the matrix Z in 
are actually the same, corresponding to the case that the two reac- 
tion networks have shared complexes (by definition only consisting of shared 
boundary chemical species; i.e., for which the corresponding columns of Z\ 
and Z2 are zero). Then one may identify the equal columns in the matrix Z, 
and thus obtain a reduced network with complex graph consisting of the union 
of the complex graphs of the two networks with the vertices corresponding to 
the shared complexes being identified. In this case the complex graph of the 
interconnected network may still be connected. 

The property of 'balancedness' is not necessarily maintained under inter- 
connection. Sufficient conditions for this to happen are given in the following 
proposition. 



Proposition 5.2. The interconnection of two balanced reaction networks 
with vectors of eqi 

if there exists (x**, x™, x™) such that 



with vectors of equilibrium constants K^^K^ 1 is again balanced if and only 



Ln 



K eq 



BJZT 






B T 2 Zl 



U 2 L b2 



Ln 



*=t= 

x 2 



(43) 



There exists such a thermodynamic equilibrium for the interconnected net- 
work if there is a partition {1, • • • ,6} = /1IJI2 such that all columns of 
BfZ^ corresponding to the index set I\ are contained in the image of Bf Zf , 
and all columns of B^Z^ corresponding to the complementary index set I2 
are contained in the image of B^Zj '. 

Proof. By assumption there exist thermodynamic equilibrium points (x J , xt^ ) 
and (xt^x^) °^ ^ ne ^ wo individual networks. That means (cf. [T4"|) ) that 



Ln(K e 1 q ) = B 1 [zT z£] Ln 



B2 [zl 



( 


x* 




* 

X bl. 


( 






* 

X b2. 



(44) 
(45) 



Now define x™ E R by taking its i-th component for i € I\ equal to the 
2-th component of x^ 2 , and for i S I2 equal to the i-th component of x\ v 
Then there exist corresponding X-^ j X 

2* such that (j4~3l) is satisfied. ■ 



Under the assumptions of Proposition 15, 2p it follows that the balanced 
interconnected network is given as 



X\ 

X-2 
X b 



-ZB 



fci(xr,*r) ° 

K 2 (x* 2 *,x* b *\ 



B T Exp I Z T Ln 



xi/x 1 
x 2 /xf 
Xb/x b _ 



(46) 



H), Z as in {i2 



'1 ! ^2 



,a;J* satisfy (03]), and /C 



where i? as in 

and K, 2 {x2* ,x b *) are proportionally related to Ki(x*,x* b ) and /C2(x 2 , x 
respectively, as follows from Proposition 13.61 



xj*) 
2b)> 



5.2 Interconnection arising from port interconnection 

The above procedure for interconnection of chemical reaction networks can 
be also interpreted as arising from power-port interconnection at the bound- 
ary chemical species. Recall the formulation (|38p of an open chemical reac- 
tion network with inputs v b being the influx/efflux of the boundary chemical 
species and Hb their chemical potentials. Then the interconnection of two 
chemical reaction networks (as above under the simplifying assumption that 
all boundary chemical species are shared) can be seen to result from the 
power-port interconnection constraints 



M61=/4>2, ^61+^62 = 0, 



(47) 



expressing that the chemical potentials of the boundary chemical species are 
equal, while the boundary fluxes of the two networks add up to zero (conser- 
vation of boundary chemical species). Indeed, the resulting interconnected 
dynamics given by the differential-algebraic equations 

-Z 1 S 1 /Ci(x*) J BfExp(Zf^ 1 (x 1 ))' 
Z 2 B 2 K 2 {x* 2 )B^ V (Z^ 2 {x 2 )) 
fJ-bi(xi) 
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x 2 _ 









+ 






-I 



[I 



Hb2(x 2 ) 



with v = Vbi = —Vb 2 acting as a vector of Lagrange multipliers, can be 
seen, after elimination of the algebraic constraint Hbi(x±) = ^,2(^2) and the 
Lagrangian multipliers v , to result in the dynamics 



x = — ZBJC(xl, x 2 )B Exp(Z n(xi,x 2 )) 



where x = (xi,X2,Xb) denotes the total vector of chemical species (internal 
species of both networks, together with the shared boundary species). Here 
the complex stoichiometric matrix Z and incidence matrix B of the intercon- 
nected network are as above, and where fc(x\,x?£) = diag(/Ci(xf ), ^(xj^)) 

and ufxi , Xo) = ^ 1 ) Xl l . Recalling that p T • 

[H2{X2)_ 

the power provided to the chemical reaction network, this implies that the 
interconnection ([171) is power-preserving, that is n^Vbi + $2^2 = 0, in line 
with the standard way of interconnecting physical networks [T9J, [181 120] • 



Recalling that fijvb is (modulo the constant RT) 



6 Model reduction of chemical reaction networks 

For many purposes one may wish to reduce the number of dynamical equa- 
tions of a complex chemical reaction network in such a way that the behavior 
of a number of key chemical species is approximated in a satisfactory way. 

Furthermore, it is important that the reduced system again allows the 
interpretation of a chemical reaction network (but with much less chemi- 
cal species involved). In fact, most of the approaches to model reduction 
of large-scale chemical reaction networks aim at doing exactly this: they 
simplify the pathways of the chemical reaction network by leaving out in- 
termediate species or complexes; see e.g. [9j[l], or the standard reduction of 
mass action enzymatic reactions to their Michaelis-Menten description |15| . 

In the following we will propose a model reduction approach which is 
based on a reduction of the complex graph. First we recall from |22| the fol- 
lowing result regarding Schur complements of weighted Laplacian matrices. 

Proposition 6.1. Consider a directed graph with vertex set V and with 
weighted Laplacian matrix L = BK{x*)B T . Then the Schur complement 
with respect to any subset of vertices V r C V is well-defined, and is the 
weighted Laplacian matrix L = BfC(x*)B T of a reduced directed graph with 
incidence matrix B whose vertex set is equal to V — V r . 

Remark 6.2. Note that the reduced graph defined by B may contain ad- 
ditional edges between the remaining vertices in V — V r . Furthermore, the 
weights of the reduced diagonal weight matrix fC(x*) are in general different 
from the weights of the original weight matrix K,(x*), and represent equiva- 
lent balanced reaction constant^j\. 



In electrical circuit theory the process of reduction of a resistive network to a resistive 
network with fewer vertices but equivalence resistance is sometimes referred to as Kron 
reduction, cf. [151 15]. 



This result can be directly applied to the complex graph, yielding a 
reduction of the chemical reaction network by reducing the number of com- 
plexes. The procedure is as follows. Consider a balanced reaction network 
described in the standard form (j2T|) 

£: x = -ZBlC(x*)B T Exp (z T Ln (Jjzj) 

Reduction will be performed by deleting certain complexes in the com- 
plex graph, resulting in a reduced complex graph with weighted Laplacian 
L = BK{x*)B T . Furthermore, leaving out the corresponding columns of 
the complex stoichiometric matrix Z one obtains a reduced complex stoi- 
chiometric matrix Z (with as many columns as the remaining number of 
complexes in the complex graph), leading to the reduced reaction network 

£ : x = -ZBlC(x*)B T Exp (Z T Ln , K(x*) > (48) 

Note that £ is again a balanced chemical reaction network governed by mass 
action kinetics, with a reduced number of complexes and with, in general, 
a different set of reactions (edges of the complex graph) . Note furthermore 
that the thermodynamic equilibrium x* of the original reaction network £ 
is a thermodynamic equilibrium of the reduced network £ as well. 

Furthermore, we can derive the following properties relating £ and £. 

Proposition 6.3. Consider the balanced reaction network £ and its reduced 
order model £ given by (f-^ffp . Denote their sets of equilibria by £, respectively 
£,. Then 8 C E . Furthermore, i/£ has deficiency zero then so does £. 

Proof. Assume that the complex graph is connected; otherwise the same 
argument can be repeated for every component (linkage class). Recall from 
(|16p that the set of equilibria £ is given as {x** | Ln T (x**)S = Ln T {x*)S}, 
where S = ZB. It follows that 8 is equivalently represented as 

£ = {x** | Ln T (x**)ZL = Ln T {x*)ZL}, 

where L := BK.(x*)B T is the weighted Laplacian matrix of the complex 
graph. Now consider a subset V r of the set of all complexes, which we wish 
to leave out in the reduced network. After permutation of the complexes we 
partition L and Z correspondingly as 

z = [z x z 2 ] 



L 



L21 L22 



where V r corresponds to the last part of the vertices. Then post-multiply 
the equations hn T (x**)ZL = Ln T \x*)ZL by the invertible matrix 



/ 



(49) 



L 22 L 2 \ I 



not changing the solution set £. It follows that 
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C {x 



ZJ [ 

Ln T (x*)ZL} = £, 



since Z = Z\ and L = L\\ — L\2L 22 L21. 

For the second statement assume that the full-order network has defi- 
ciency zero, i.e., kei Z f^\imB = 0, or equivalently kerZP|imL = 0. Post- 
multiplication of L with the matrix in (149)) yields 



which implies kerZP|imL = 0, i.e., zero-deficiency of the reduced-order 
network E. ■ 

The precise properties of this model reduction approach are currently 
under investigation. At least the approach seems to have favorable prop- 
erties for the shortening of linear pathways, leaving out the intermediary 
complexes. Furthermore, in most cases it seems natural not to leave out 
complexes containing boundary chemical species. Note that the model re- 
duction approach results in the reduction of a number of chemical species in 
the reduced network once these chemical species are only contained in the 
complexes that are deleted in the reduction procedure. 

6.1 Reversible enzymatic reaction 

As a very simple example, let us apply the above procedure to a standard 
reversible enzymatic reaction represented as 




= 



E + X 



where X, Y denote chemical species, E the enzyme, and / the intermediate 
species (or complex). The three complexes involved are E + X, I, E + Y 
and we want to leave out the intermediate complex /. Since the transposed 
stoichiometric matrix S T has full row-rank the Wegscheider conditions (|15p 
are trivially satisfied, guaranteeing the existence of a thermodynamic equi- 
librium. Denote for later use the concentrations of the species X, Y, E by 
x,y,e, and their thermodynamic equilibrium concentrations by x*,y*,e*. 

Furthermore, denote the balanced reaction constant of the first reversible 
reaction by kx , and of the second one by Ky . Then the weighted Laplacian 
matrix of the complex graph is given by 

k x -k x 

L = —Kx Kx + Ky —Ky 
— Ky Ky 

Deletion of the intermediate complex / corresponds to taking the Schur 
complement of L with respect to the second row and column, yielding the 
reduced weighted Laplacian matrix 

L = 

where the reduced balanced reaction constant k is computed as 

K X Ky 
KX + Ky 

Furthermore, the reduced complex stoichiometric matrix Z expressing the 
remaining complexes E + X, E + Y in terms of the species X, Y, E is given 

bB 

"1 0" 



z 



It follows that the reduced dynamics P5|) is given as 
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18 Strictly speaking, the matrix Z should be complemented by a zero row correspond- 
ing to the chemical species 7, which is however not present anymore in the remaining 
complexes. 



which is in mass action kinetics form, and has (x*,y*,e*) as a thermo- 
dynamic equilibrium. Note that in the reduced model the total enzyme 
concentration e is preserved (as to be expected). 

It is easily checked that the equations for the x, y and e values of the 
thermodynamic equilibria of the reduced network are the same as for the 
original network. On the other hand, in contrast with the thermodynamic 
equilibrium equations for i of the original model, the reduced model allows 
any value of i as equilibrium value; in accordance with Proposition 16.31 
Furthermore, both networks are trivially zero-deficient. 

7 Conclusions 

In this paper we have provided a compact, geometric, formulation of the 
dynamics of mass action chemical reaction networks possessing a thermo- 
dynamic equilibrium. This formulation clearly exhibits both the structure 
of the complex graph and the stoichiometry, and furthermore admits a di- 
rect thermodynamical interpretation. Exploiting this formulation we were 
able to recover, for this subclass of mass action kinetics chemical reaction 
networks, some of the results in the fundamental work \\.4\ [T3l [TOl [TT] in a 
simple and insightful way, without having to rely on the properties of defi- 
ciency zero or one. Drawing inspiration from |191 [18, 20J, we have shown how 
the framework leads to an input-state-output formulation of open chemical 
reaction networks, and how this may be used for interconnection. Further- 
more, we have indicated how this formulation, in particular its emphasis on 
the Laplacian matrix of the complex graph, may be used for a systematic 
model reduction procedure based on the elimination of certain intermediate 
complexes. 

Current research is taking place in a number of directions. We believe the 
striking similarity with consensus dynamics can be further exploited. Also, 
the use of our formulation for the analysis of steady states corresponding to 
non-zero (constant) boundary fluxes is under study, taking into account the 
possibility of multiple steady states [6] ; see also the work [3] on input-output 
stability. The ramifications of the new model reduction procedure are being 
investigated, including the comparison with classical reduction schemes like 
(reversible) Michaelis-Menten reaction rates; see e.g. [15]. Furthermore, the 
Laplacian matrix can not only be used for model reduction, but also for 
other analysis purposes of the complex graph, e.g. for its decompostion pQ. 
Finally, a challenging question is the extension to biochemical reaction net- 
works, where the reactions are mostly enzymatic reactions, and the amount 



of enzymes, through the activity of the gene and protein networks, will 
depend on e.g. the concentration of a number of chemical species (metabo- 
lites). This will lead to the introduction of an additional network structure 
originating from the regulatory, and possibly competing, feedback loops. 
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